﻿#include <iostream>
#include <opencv2\opencv.hpp>

//https://blog.csdn.net/qq_34745295/article/details/79411564

//4. 使用米粒图像，分割得到各米粒，
//首先计算各区域(米粒)的面积、长度等信息，
//进一步计算面积、长度的均值及方差，分析落在3sigma范围内米粒的数量。

using namespace cv;
using namespace std;


void main()
{
	Mat img_src, gray, erosion, dilation, backImg, rice, img_otsu, img_seg;
	
	//cout << element << endl;

	const char *fn = "rice.jpg";
	img_src = imread(fn);

	if (img_src.empty())
	{
		cout << "无法打开源图像." << endl;

		system("pause");
		return ;
	}
	//resize(img_src, img_src, cv::Size(img_src.rows*2, img_src.cols*2));//用原来大小和X2倍后检测出来的的米粒数和3sigma合格数比例都不太一样
	cvtColor(img_src, gray, COLOR_BGR2GRAY);

	//形态学处理，去除噪声	  构造模板，5次腐蚀  开运算：先腐蚀后膨胀
	//morphologyEx(bw, bw, MORPH_OPEN, element);
	//morphologyEx(bw, bw, MORPH_CLOSE, element);
	Mat kernel = getStructuringElement(MORPH_RECT, Size(7, 7)); //构造卷积核
	erode(gray, erosion, kernel, Point(-1, -1), 5); //进行5次腐蚀
	dilate(erosion, dilation, kernel, Point(-1, -1), 5);//进行5次膨胀

	//原图(灰度图)减去背景得到米粒形状
	backImg = dilation;
	rice = gray - backImg;
	imshow("背景图", backImg);
	imshow("灰度图-背景图", rice);

	// 大津算法阈值化
	threshold(rice, img_otsu, 0, 0xff, THRESH_OTSU);
	imshow("背景图+OTSU", img_otsu);


	// 以下是图像分割
	img_seg = img_otsu.clone();
	vector<vector<Point>> cnts;
	findContours(img_seg, cnts, RETR_LIST, CHAIN_APPROX_SIMPLE);
	//**#轮廓检测，注意输出结果应有3个，其他著作上为2个输出结果，亲测错误**
	//ret1, contours, hierarchy = cv2.findContours(ret1, cv2.RETR_LIST, cv2.CHAIN_APPROX_SIMPLE)


	//统计变量
	int count = 0;          //米粒数量
	//面积变量
	float area_sum = 0.0;   //米粒面积求和
	float v_EX_area = 0.0;  //米粒面积平均
	float v_DX_area = 0.0;  //米粒面积方差
	float v_Standard_DX_area = 0.0;  //米粒面积标准差
	float v_area2_sum = 0.0;//米粒求和（每个米粒面积的平方）
	//长度变量
	float length_sum = 0.0;   //米粒长度的求和
	float v_EX_length = 0.0;  //米粒长度的平均值，期望，作为mu
	float v_DX_length = 0.0;  //米粒长度的方差，作为sigma
	float v_Standard_DX_length = 0.0;  //米粒长度标准差
	float v_length2_sum = 0.0;//米粒求和（每个米粒长度的平方）
	//
	float check_3sigma_length[2];//3sigma的范围(米粒长度)
	float check_3sigma_area[2];//3sigma的范围(米粒面积)
	//缓存每个米粒长度和面积
	list<float> list_lengths;
	list<float> list_areas;
	//符合3sigma的米粒数量
	int legal_3sigma_length_count = 0;
	int legal_3sigma_area_count = 0;

	for (int i = cnts.size() - 1; i >= 0; i--)
	{
		int length = 0;
		vector<Point> cnt = cnts[i];
		float area = contourArea(cnt);//轮廓面积
		double dis_max = 0.0;//米粒长度.  这里枚举cnt所有点的距离取最大.    用矩阵包围盒算？
		Point dis_max_point0 = cnt[0];
		Point dis_max_point1 = cnt[0];

		//计算米粒长度
		for (vector<Point>::iterator it1 = cnt.begin(); it1 != cnt.end(); it1++)
		{
			Point point1 = *it1;
			for (vector<Point>::iterator it2 = cnt.begin(); it2 != cnt.end(); it2++)
			{
				Point point2 = *it2;
				if (point1.x == point2.x && point1.y == point2.y)
				{

				}
				else
				{
					double distance;
					distance = powf(((float)point1.x - (float)point2.x), 2) + powf(((float)point1.y - (float)point2.y), 2);
					distance = sqrtf(distance);
					if (distance >= dis_max)
					{
						dis_max = distance;
						dis_max_point0 = point1;
						dis_max_point1 = point2;
					}
				}
			}
		}


		if (area < 10) // 滤除面积小于10的分割结果：可能是噪声
		{
			cout << "idx(过滤噪点):" << i << "  米粒面积: " << area << " 米粒长度:" << dis_max << endl;
			vector<vector<Point>> tmp_cnts;
			tmp_cnts.push_back(cnt);
			drawContours(img_src, tmp_cnts, -1, Scalar(0, 0, 255), 1);

			line(img_src, dis_max_point0, dis_max_point1, Scalar(255, 0, 0), 1);

			putText(img_src, to_string(i), Point((dis_max_point1.x - dis_max_point0.x) / 2 + dis_max_point0.x, (dis_max_point1.y - dis_max_point0.y) / 2 + dis_max_point0.y), CV_FONT_HERSHEY_PLAIN, 2.5, Scalar(0, 0xff, 0xff));

	

		}
		else
		{
			count++; // 统计米粒数量
			area_sum += area;

			list_lengths.push_back(dis_max);
			list_areas.push_back(area);

			cout << "idx:" << i << "  米粒面积: " << area << " 米粒长度:" << dis_max << endl;
			length_sum += dis_max;
			v_length2_sum += (dis_max * dis_max);

			v_area2_sum += (area * area);

			vector<vector<Point>> tmp_cnts;
			tmp_cnts.push_back(cnt);
			drawContours(img_src, tmp_cnts, -1, Scalar(0, 0, 255), 1);
			line(img_src, dis_max_point0, dis_max_point1, Scalar(255, 0, 0), 1);
			putText(img_src, to_string(i), Point((dis_max_point1.x - dis_max_point0.x) / 2 + dis_max_point0.x, (dis_max_point1.y - dis_max_point0.y) / 2 + dis_max_point0.y), CV_FONT_HERSHEY_PLAIN, 0.5, Scalar(0, 0xff, 0));
		}
	}

	v_EX_length = length_sum / count;
	v_DX_length = v_length2_sum / count - v_EX_length * v_EX_length;
	v_Standard_DX_length = sqrtf(v_DX_length);

	v_EX_area = area_sum / count;
	v_DX_area = v_area2_sum / count - v_EX_area * v_EX_area;
	v_Standard_DX_area = sqrtf(v_DX_area);


	cout << ">10面积的米粒数量： " << count <<  endl;
	cout << "米粒面积均值:" << v_EX_area << "   米粒面积方差:" << v_DX_area << "  米粒面积标准差:" << v_Standard_DX_area <<  endl;
	cout << "米粒长度均值:" << v_EX_length << "   米粒长度方差:" << v_DX_length << "  米粒长度标准差:" << v_Standard_DX_length << endl;

	//对于米粒长度，横轴取值u-3sigama到u3sigama之间的范围
	check_3sigma_length[0] = v_EX_length - 3 * v_Standard_DX_length;
	check_3sigma_length[1] = v_EX_length + 3 * v_Standard_DX_length;
	//对于米粒面积，横轴取值u-3sigama到u3sigama之间的范围
	check_3sigma_area[0] = v_EX_area - 3 * v_Standard_DX_area;
	check_3sigma_area[1] = v_EX_area + 3 * v_Standard_DX_area;


	for (list<float>::iterator it = list_lengths.begin(); it != list_lengths.end(); it++)
	{
		float length = *it;
		if (length >= check_3sigma_length[0] && length <= check_3sigma_length[1])
		{
			legal_3sigma_length_count++;
		}
	}
	for (list<float>::iterator it = list_areas.begin(); it != list_areas.end(); it++)
	{
		float area = *it;
		if (area >= check_3sigma_area[0] && area <= check_3sigma_area[1])
		{
			legal_3sigma_area_count++;
		}
	}
	cout << "3sigma范围(长度)： " << check_3sigma_length[0] << "  " << check_3sigma_length[1] << "  合格数:" << legal_3sigma_length_count << endl;
	cout << "3sigma范围(面积)： " << check_3sigma_area[0] << "  " << check_3sigma_area[1] << "  合格数:" << legal_3sigma_area_count << endl;

	imshow("原图", img_src);

	waitKey(0);
}

